Micro-structure of damage in thermally activated fracture of Lennard-Jones systems 
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We investigate the effect of thermal fluctuations on the critical stress and the micro-structure 
of damage preceding macroscopic fracture of Lennard-Jones solids under a constant external load. 
Based on molecular dynamics simulations of notched specimens at finite temperature, we show 
that the crystalline structure gets distorted ahead of the crack in the secondary creep regime. The 
damage profile characterizing the spatial distribution of lattice distortions is well described by an 
exponential form. The characteristic length of the exponential form provides the scale of damage 
which is found to be an increasing function of the temperature: At low temperature damage is 
strongly localized to the crack tip, while at high temperature damage extends to a broader range 
leading to more efficient relaxation of overloads. As a consequence, the stress intensity factor 
decreases with increasing temperature. The final macroscopic failure of the system occurs suddenly 
which is initiated by the creation of vacancies and voids. The creep strength exhibits the inverse 
square root scaling with the notch size corrected by the extension of the process zone. 

PACS numbers: 62.20.M-,46.50.+a,81.40.Np 



I. INTRODUCTION 

Understanding sub-critical fracture occurring under 
constant or periodic external loads below the fracture 
strength of materials is a very important scientific prob- 
lem with a broad spectrum of technological applica- 
tions. Depending on materials' details several micro- 
scopic mechanisms may contribute to the time dependent 
response and final failure of samples. In spite of this di- 
versity, experiments and theoretical investigations have 
revealed that sub-critical fracture obeys scaling relations 
both on the macro and micro scales p], Q ■ 

Thermal activation of micro-crack nucleation and 
growth is one of the primary mechanisms underlying sub- 
critical fracture. Theoretical and experimental investi- 
gations have shown that both in homogeneous and het- 
erogeneous materials thermally driven stress fluctuations 
can be responsible for the finite lifetime of loaded speci- 
mens [3, 4] leading to the emergence of universal scaling 
laws such as the Andrade law and time-to-failure power 
laws [a, [6[, furthermore, to the scalin g be havior of wait- 
ing times between micro- fractures |7H11|. Experiments 
have also revealed that thermally activated slow crack 
advancements affect even the surface roughness of grow- 
ing cracks [12j . 

In order to understand the mechanisms of the sub- 
critical fracture of crystalline and glassy materials at the 
atomistic scale, molecular dynamics simulations [l3j per- 
formed at finite temperatures proved to be indispensable 
[l3 - [l9| . Computer simulations of Lennard-Jones parti- 
cle systems made possible to understand the origin of the 
Griffith-type crack nucleation at the atomic scale show- 
ing the im por tance of the creation of vacancy clusters 
and voids [15l - l22 |. The simulation technique was also 
extended to study dynamic fracture phenomena (23| . 

In the present paper we study thermally activated sub- 
critical fracture using molecular dynamics simulations of 
a Lennard-Jones particle system in two dimensions. A 



rectangular sample with a triangular lattice structure is 
constructed which is then subject to a constant external 
load at different notch length varying also the tempera- 
ture. We show that the critical stress of the system has 
an inverse square root dependence on the notch length as 
expected from linear fracture mechanics (LFM) corrected 
by the extension of the fracture process zone. The stress 
intensity factor at a given notch length decreases with 
increasing temperature due to the more efficient relax- 
ation of overloads in the vicinity of the crack. In order 
to characterize the micro-structure of the loaded sam- 
ple, we introduce the six-fold order parameter and an- 
alyze its spatial distribution during the secondary creep 
regime. Computer simulations revealed that damage con- 
centrates ahead of the crack tip and decays to a nearly 
homogeneous background at larger distances. The dam- 
age profile has an exponential form from which the char- 
acteristic scale of damage can be determined. Increasing 
the temperature the damage profile broadens and the 
background damage increases which are consistent with 
the decreasing stress intensity factor. The results are dis- 
cussed in comparison to the micro-structure of heteroge- 
neous solids with only quenched disorder under quasi- 
static loading just before ultimate failure [24j . 



II. MODEL 

We perform molecular dynamics (MD) simulations of 
two dimensional solids which consists of particles. The 
particles interact through the Lennard-Jones (LJ) poten- 
tial 



U{r) = 4e 



(1) 



where r, e and a are the distance between particles, the 
strength of the potential and the characteristic length 
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FIG. 1. (Color online) Schematic picture of the system: The 
sample is loaded through springs in the vertical direction. 
Bottom and upper sides are attached to the Nose-Hoover ther- 
mostat to control the temperature. In the horizontal direction 
periodic boundary condition is applied. 



FIG. 2. Time evolutions of a sample at the temperature T = 
0.1 and the load a = 1.45 with the initial crack length 2L = 
20: (a) an initial state, and (b) a configuration after uploading 
and equilibration, (c) As time evolves, the lattice structure 
gets distorted, the crack grows and voids and defects nucleate 
ahead of the crack tip. (d) The material's fracture occurs 
suddenly following a long evolution process during which the 
overall deformation hardly changes. 



of the particle, respectively. For the computational ef- 
ficiency, we cut the interaction potential Eq. (fTJ) at the 
distance 3a when calculating the inter-particle force, i.e. 
if the distance r between two particles exceeds 3a, the in- 
teraction disappears. The parameters e, a and the mass 
of particles m are set to be unity in all the simulations. 
The temperature T is defined as the average kinetic en- 
ergy of the system T = Yli=i mv?/2NkB, where de- 
notes the velocity of i-th particle, N is the number of 
particles, and fee is the Boltzmann's constant. Hereafter 
we also set the Boltzmann's constant fee to be unity. 

A schematic picture of the computational geometry is 
presented in Fig. [TJ In order to control the tempera- 
ture, Nose- Hoover thermostats are attached to two parti- 
cle layers on the top and bottom of the system [lEHBl- In 
the horizontal direction we apply periodic boundary con- 
dition to make the crack advance from the initial crack in 
the middle of the system and to prevent the vaporization 
of particles due to the thermal effect. 

The initial state of the sample is prepared as follows: 
The particles are placed on a regular triangular lattice at 
zero temperature T = then the temperature of the sys- 
tem is slowly increased up to the desired value. We take 
the width of the system as the corresponding equilibrium 
value for the temperature. Simulations were carried out 
with a fixed value of the width Iq ~ 100 of the lattice. 
As a crucial component of the model, an initial crack is 
created by removing four particle layers with length 2L 
in the middle of the system. In this way the crack has a 



nearly rectangular shape with extensions 2L x 3.5a. Note 
that the distance 3.5a of the two sides is larger than the 
cutoff length 3a of the potential, which ensures that the 
removed particles form a crack in the model. In order to 
prevent the crack from healing due to thermal fluctua- 
tion, the particles on the upper and bottom crack faces 
interact through the repulsive part of the LJ potential 
Eq. ([T]). After equilibrating the system, the load is ap- 
plied in such a way that the bottom and top particle 
layers are attached to loading bars through spring ele- 
ments. The bottom bar is fixed while the upper one can 
move vertically controlling the load a on the system. The 
bars are treated as rigid bodies with a mass 10 -3 per unit 
length, while the loading springs have zero mass and a 
spring constant 10. The stiffness of the solid is about 
83 for the parameters used in the simulations, which is 
significantly larger than the one of springs. 

In order to suppress the disturbing effect of elastic 
waves generated by the external loading process, the 
loading bar moves from the initial position until the de- 
sired stress value a is reached at a relatively low increas- 
ing rate da/dt — 10~ 2 . Simulations showed that the 
value of the loading rate slower than da/dt = 10~ 2 did 
not affect the time evolution of the system. 
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FIG. 3. (Color online) Vertical size l(t) of the system as a 
function of time for several temperatures with the initial crack 
length 2L — 20. Elastic waves generated by the loading pro- 
cess and the void creations due to the thermal noise cause the 
small oscillations of l(t). The vertical arrows indicate the time 
window for the sample of T = 0.2 where the damage parame- 
ter d is evaluated (see Sec.|V}. Applied loads are taken to be 
a = 2.05, 1.7, 1.45, 1.35, and a = 1.1 for T = 0, 0.06, 0.1, 0.14, 
and T = 0.2, respectively. 



III. TIME EVOLUTION 

The time evolution of the solid subject to a constant 
load is calculated numerically by the second-order sym- 
plectic integrator method with the time step dt = 10~ 3 , 
while the time evolution of the thermostat and the bar is 
obtained using the Euler method. Figure[2]presents snap- 
shots of the time evolution of the system at the temper- 
ature T = 0.1, the load a — 1.45, the initial crack length 
2L = 20 and the width of the system Iq ~ 100. Starting 
from the initial state (Fig.^a)) we gradually increase the 
external load a up to the desired value and equilibrate the 
system to achieve the desired temperature. Figure &b) 
shows the sample after the ramp loading, i.e., when the 
stress becomes constant for the rest of the simulation. It 
is observed in Fig. [^c) that the crystalline structure of 
the lattice gets distorted. Nucleation of defects and voids 
due to thermal fluctuations is also observed, which facil- 
itates the advancement of the crack. When the specimen 
is separated into two parts by a crack like in Fig. [U[d) , 
the system is considered to be fractured. 

On the macroscopic level, the time evolution of the 
system can be monitored by measuring the vertical size 
l(t) of the sample in the direction of the external load 
as a function of time t. This quantity is presented in 
Fig. [3] for several temperatures with the initial crack 
length 2L = 20. Note that the sample of Fig. [2] cor- 
responds to the curve of T = 0.1 in Fig. [3] The plateau 
regime of the deformation is the consequence of the slow 
dynamics of the system, where the lattice structure gets 
gradually distorted while the length of the crack hardly 



changes. The final macroscopic failure occurs rapidly, 
which is characterized by the sudden acceleration of the 
vertical size l(t) and by its diverging derivative in Fig. |3] 
The acceleration of the vertical size is preceded by the 
creation of defects and voids due to thermal activation. 
Small oscillations are observed along the plateau of l(t), 
which are caused by elastic waves generated by the load- 
ing process and the void creations due to the thermal 
noise. 



IV. MACROSCOPIC STRENGTH 

The main difficulty of the investigation of thermally 
activated fracture of LJ systems is the enormous com- 
putational time required due to the slow dynamics of 
the system. To overcome this problem, we introduce 
a threshold time of the simulation; After the stress a 
reaches the desired value, we limit the simulation time 
to the finite threshold time t t h- This threshold gives a 
lower limit of the lifetime of the sample; If no complete 
failure is achieved by following the time evolution of the 
system up to t t h, the lifetime of the sample is assumed 
to be infinite at the given external load a. The highest 
load at which no macroscopic failure occurs before tth 
defines the critical load a c of the system which separates 
the regimes of finite and infinite lifetimes. The critical 
load a c of the system is determined by a sequence of sim- 
ulations performed with a constant load increment 0.05, 
which provides sufficient accuracy for a c . In this paper, 
the threshold time is fixed to be tth = 10 3 . 

In order to investigate the effect of the notch size 2L 
and of the temperature on the macroscopic strength er c , 
we carry out simulations with the following parameters: 
2L = 0,5,10,15,20,25,30,35,40, and 2L = 45 for the 
temperatures T = 0,0.06,0.1,0.14, and T = 0.2. We 
do not show the data of 2L = 0, 5 at T — because 
the crack did not grow at the crack tip. The value of 
<7 C is determined by averaging over 10 samples at each 
parameter set. 

Based on linear fracture mechanics we expect the in- 
verse square root dependence of the critical load a c on 
the notch size L as it has been proved by an analytical 
calculation [27j. In order to capture the disordering ef- 
fect of temperature, we introduce a characteristic length 
£ which modifies the notch size into an effective one L + £. 
Then we write a c in the following form 

a (L t T) = -^ + c (2) 

where K denotes the stress intensity factor [24| H3| and 
c is a finite size correction of the system width. The 
finite size correction c is a negative value and its system 
size dependence is plotted in the inset of Fig. S) It is 
important to emphasize that the value of the correction 
factor c goes to zero with increasing system size Iq. To 
numerically demonstrate the validity of the functional 
form Eq. ©, in Fig.0]we plot (1 / a' c (L)) 2 , where cr' c (L) = 
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FIG. 4. (Color online) Crack length dependence of the crit- 
ical stress: (l/a' c (L)) , where cr' c (L) = a c (L) — c and c is a 
finite size correction, is plotted as a function of L for several 
temperatures averaging over 10 samples. The error bars are 
smaller than the symbol size. The dotted lines are obtained 
by fitting straight lines in a good agreement with Eq. ©• De- 
viations near L > 20 of the data from the linear fit occur due 
to finite size effects, while the deviations for the small notch 
size indicate the dominance of thermal noise in the fracture 
process. Inset: System width lo dependence of the finite size 
correction c. 



cr c (L) — c, as a function of half of the notch length L. It is 
observed in Fig. @]that straight lines are obtained with a 
good accuracy at all temperatures by fitting. Deviations 
from the linear form are observed for both large and small 
values of the initial crack length: In the large crack length 
regime 20 < L, deviations occur due to the finite size 
of the system, on the other hand, in the limit of small 
cracks L < 2.5, the stress concentration at the crack tip 
becomes less dominating so that the fracture process is 
mainly controlled by thermal noise emerging in the entire 
volume of the system. By fitting these data from L = 2.5 
to 20, we estimate the stress intensity factor K and the 
characteristic length £ of Eq. (T2|) for each temperature. 

Figs.[5Ja) and (b) present the temperature dependence 
of the stress intensity factor K and of £, respectively. It 
is observed that the stress intensity factor K decreases 
with increasing temperature T, while the value of £ is an 
increasing function of T. In Fig.[SIa), at higher tempera- 
ture thermally activated particle motion has a higher in- 
tensity which facilitates the relaxation of the system and 
gives rise to a decreasing K. In Fig. [5jb) the results im- 
ply that at higher temperature a broader zone is formed 
ahead of the crack tip where distortion of the crystalline 
structure and plastic deformation can occur due to the 
interplay between the stress field and the thermal noise. 
The length scale £ can be interpreted as the linear exten- 
sion of the fracture process zone (FPZ) ahead of the crack 
tip which appears due to the disordering effect of ther- 
mal noise on the lattice structure. In the vicinity of the 
crack tip the high stress concentration makes the system 




FIG. 5. (a) Temperature dependence of the stress intensity 
factor K: As the temperature increases the value of K is 
monotonically decreasing, (b) Solid line shows the temper- 
ature dependence of the characteristic length £ defined by 
Eq. (fj?| . Dashed line represents the extension of the distorted 
zone £fpz determined by the analysis of the micro-structure 
of the system (This analysis is carried out in Sec.[V]). Simu- 
lations are carried out with a fixed notch size 2L = 20. 




FIG. 6. (Color online) Spatial distribution of the bond- 
orientational order parameter averaged on a rectangular 2x2 
mesh: The temperature and the initial crack length are taken 
to be T — 0.06 and 2L — 20. For the sake of clarity, we present 
a magnified view of the vicinity of the crack. The color code 
indicates the modification of the triangular-crystalline struc- 
ture due to the high stress concentration and to the thermal 
noise. Deep blue indicates the order parameter less than 0.9. 
In this figure, the blue region corresponds to the crack. 

more sensitive to thermal fluctuations. As a consequence, 
a distorted zone emerges which affects the scaling of <7 C 
with the notch size L. Similar effect is observed in het- 
erogeneous quasi-brittle materials where damage in the 
form of micro-cracks is concentrated ahead of the crack 




V. MICRO-STRUCTURE OF DAMAGE 

In order to obtain a direct quantitative measure of 
the extension of the process zone formed ahead of the 
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crack tip, we carry out a detailed analysis of the micro- 
structure of the system before the acceleration of defor- 
mation starts. To characterize the crystalline order, we 
introduce the bond-orientational order parameter (f>k de- 
fined for particle k as follows [28[ 



(3) 



Here nk denotes the number of neighbors of particle k and 
Oki is an angle between a fixed axis such as the x axis and 
the bond connecting particle k to particle I. The value 
of ipk can vary between zero and one such that <pk = 
1 quantifies a perfectly ordered triangular lattice, while 
(j)k ~ is obtained for a disordered fluid-like system. 
Intermediate values of the order parameter characterize 
the degree of distortion of the lattice structure. 

In order to characterize the spatial variation of the 
micro-structure, we calculate a locally averaged bond- 
orientational order parameter on a coarse-grained mesh 
in the sample. The mesh size is taken to be 2 x 2. The 
order parameter is evaluated for a single sample, but 
30 snapshots are taken from the plateau regime of the 
deformation-time diagrams in order to improve the statis- 
tics. As shown in Fig. [3l the two vertical arrows indicate 
the time window of the plateau regime for the case of 
temperature T — 0.2. A typical example of the order 
parameter <p for the temperature T = 0.06 is presented 
in Fig. [6] It is observed that the crystalline structure has 
strong distortion in the vicinity of the crack tip where 
the stress concentration is high. Away from the crack 
tip, and even near the middle of the crack, the order pa- 
rameter is homogeneously distributed and it has a large 
value close to unity, i.e. the crystalline structure is re- 
tained there. 

To obtain a detailed quantitative measure of the dis- 
ordering effect of the stress concentration assisted by the 
thermal noise, we introduce a damage parameter with 
the following definition 



4 = 1- 



(4) 



which has a large value dk ~ 1 at locations where the 
lattice structure is most distorted. Based on the above 
meshing technique, we average dk on each plaquette of 
the mesh, then make a projection on the horizontal axis, 
and determine the damage parameter along the initial 
crack axis as a function of the distance from the crack 
tip x; At each position x, the averaged damage param- 
eter is averaged again on the vertical column of 5 pla- 
quettes. The damage profile d(x) obtained by this way 
is presented in Fig. [7] for several different temperatures. 
Note that the calculations are performed at the same load 
value a = 1.0 for all the temperatures at the given notch 
length 2L = 20. It is observed that d(x) monotonically 
decreases with increasing distance x from the crack tip 
as it is expected, and converges to finite values at long 
distances. The figure presents that the profiles can be 
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FIG. 7. (Color online) Damage profile, i.e. the damage pa- 
rameter d as the function of the distance x from the crack tip 
for the initial crack length 2L = 20 at several temperatures. 
The damage profile exhibits an exponential decay with the 
form Eq. ([5| to a uniform background damage. All the three 
parameters characterizing the functional form of d(x) depend 
on the temperature (see Fig. 0(b) and Fig. [8]). 




FIG. 8. The parameters A and B of the damage profile as 
function of temperature for the notch size 2L = 20. The 
value of the background disorder A rapidly increases with the 
temperature, while the multiplication factor B has a moderate 



well described by the exponential form 

d(x) = A + Bexp (—x/^fpz) 



(5) 



where the parameter A represents the uniform back- 
ground damage, £,fpz is the characteristic length of the 
extension of the distorted zone, and the multiplication 
factor B represents the maximum value of the additional 
damage above the background achieved at the crack tip. 
The functional form Eq. ([5]) provides an excellent fit of 
the numerical data in all cases. Therefore we determine 
the temperature dependence of the parameters A, B, and 
£,fpz (see Fig. [5] and Fig. EJb)). Note that the dam- 
age parameter d at zero temperature decreases to zero at 
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large distances, while it takes non-zero value at the crack 
tip. This clearly indicates that the distortion of the lat- 
tice is caused by the high stress concentration arising in 
the vicinity of the crack tip. Increasing the temperature, 
the thermal noise gets stronger which lightens the stress 
concentration near the crack tip by distorting the lat- 
tice structure. This mechanism has the consequence that 
with increasing temperature the localization of stress at 
the crack tip decreases which is quantified by the increas- 
ing characteristic length £_fpz ( se e Fig.[5jb) and Fig. 

The most important characteristics of the micro- 
structure is represented by the length scale £fpz which 
provides a measure of the extension of the highly dis- 
torted zone. It can be observed in Fig. [5jb) that ^fpz 
increases when the temperature gets higher. It shows 
that increasing thermal noise in the system gives rise to 
a larger extension of the fracture process zone. It can 
also be observed in Fig. [5jb) that £ and £p pz have a lin- 
ear relation which implies that the micro-structure based 
length £p pz provides a good measure of the macroscopic 
fracture process zone. 

VI. CONCLUSIONS AND DISCUSSION 

We carried out a detailed study of the creep rupture 
of Lennard- Jones particle systems by means of molecular 
dynamics simulations. In order to control the tempera- 
ture a Nose-Hoover thermostat was coupled to the sample 
which was then subject to a constant load through spring 
elements. The main focus of our work was to understand 
the effect of thermal fluctuations on the creep strength 
and to investigate the micro-structure of the evolving sys- 
tem before the onset of rapid crack growth leading even- 
tually to macroscopic failure. Simulations showed that 
the critical stress of samples has an inverse square root 
dependence on the notch length corrected by the exten- 
sion of the fracture process zone. The fracture process 
zone arises due to the distortion of the lattice structure 
as a consequence of the interplay of the high stress con- 
centration at the crack tip and of the thermal noise. The 
final failure occurs rapidly after the creation of vacancies 
and voids in the process zone. 

In order to obtain a direct quantitative measure of the 
extension of the fracture process zone, we introduced a 
damage parameter and analyzed its spatial distribution. 
Our calculations revealed that damage concentrates in 
the vicinity of the crack tip and decays exponentially to 
a homogeneous background damage level. In spite of the 
inverse square root form of the stress distribution ahead 



of the crack expected from linear fracture mechanics, the 
damage profile is described by an exponential form which 
let us introduce a characteristic length scale of damage 
£fpz- The highly distorted regime of extension £fpz is 
analogous to the fracture process zone observed in various 
types of systems. 

A similar exponential form of the damage profile was 
recently found for heterogeneous materials where ther- 
mal activation does not play a role [24j. The authors 
investigated quasi-brittle fracture processes in the frame- 
work of the fuse model gradually increasing the load on 
notched samples up to failure. Evaluating the density of 
micro-cracks in the critical state of the sample just before 
macroscopic failure, it was found that damage concen- 
trates ahead of the crack. The damage profile proved to 
have an exponential form similar to our one. It was ar- 
gued that micro-cracks concentrate at the crack tip due 
to the high stress concentration, however, as a counter 
effect they screen the stress field which then leads to the 
exponential decay instead of the inverse square root form 

A very interesting question addressed by Ref. [24j is 
the crossover between the two regimes where the fracture 
strength of the system is dominated by stress concentra- 
tion at the notch and by the disorder, respectively. The 
deviation from the linear behavior in Fig. U at the limit 
of small notch sizes already indicates the dominance of 
thermal noise in the fracture process. However, due to 
technical reasons we could not further decrease the ini- 
tial crack length which prevents us to deduce quantitative 
conclusion on the crossover phenomenon. 
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